#### installing packages if not installed ####

list.of.packages <- c('tidyverse',
                      'ggplot2',
                      'lubridate',
                      'estimatr',
                      'ggthemes',
                      'readr',
                      'readxl',
                      'rdrobust',
                      'gridExtra')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#### attaching libraries ####

suppressPackageStartupMessages(
  
  {
    
    library(tidyverse)
    library(meta)
    library(ggplot2)
    library(lubridate)
    library(estimatr)    
    library(ggthemes)
    library(readr)
    library(readxl)
    library(rdrobust)
    library(texreg)
    library(gridExtra)
    
  }
  
)

#### loading in and generating hit rate data #### 

# Austin
# 1) Traffic Stop Hit Rates x

# Los Angeles
# 1) Pedestrian Hit Rates x
# 2) Traffic Hit Rates x

# Philly 
# 1) Pedestrian Hit Rate
# 2) Traffic Hit Rate

# Seattle 
# 1) Terry Stop Hit Rate x

# Austin Traffic Stop Data

load(file = "data/rpf.RData")

hit_rate_traffic_aus = rpf %>% 
  mutate(count = 1) %>% 
  mutate(hit = search_found,
         hit_alc = search_found_alc) %>% 
  mutate(arrest_sum = arrest) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            search_found_cash_m = mean(search_found_cash, na.rm = TRUE),
            search_found_substances_m = mean(search_found_substances, na.rm = TRUE),
            search_found_weapons_m = mean(search_found_weapons, na.rm = TRUE),
            search_found_cash_s = sum(search_found_cash, na.rm = TRUE),
            search_found_substances_s = sum(search_found_substances, na.rm = TRUE),
            search_found_weapons_s = sum(search_found_weapons, na.rm = TRUE),
            hit = mean(search_found, na.rm = TRUE),
            hit_alc = mean(hit_alc, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest_sum, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

hit_rate_traffic_aus$blmrv = 
  hit_rate_traffic_aus$date - 
  min(hit_rate_traffic_aus$date[hit_rate_traffic_aus$blm == 1])

hit_rate_traffic_aus$monthfe = substring(hit_rate_traffic_aus$date, 6, 7)

hit_rate_traffic_aus$search_found_cash_m = 
  hit_rate_traffic_aus$search_found_cash_m / 
  max(hit_rate_traffic_aus$search_found_cash_m, na.rm = TRUE)
hit_rate_traffic_aus$search_found_weapons_m = 
  hit_rate_traffic_aus$search_found_weapons_m /
  max(hit_rate_traffic_aus$search_found_weapons_m, na.rm = TRUE)
hit_rate_traffic_aus$search_found_substances_m = 
  hit_rate_traffic_aus$search_found_substances_m /
  max(hit_rate_traffic_aus$search_found_substances_m, na.rm = TRUE)

hit_rate_traffic_aus$search_found_cash_s = 
  hit_rate_traffic_aus$search_found_cash_s / 
  max(hit_rate_traffic_aus$search_found_cash_s, na.rm = TRUE)
hit_rate_traffic_aus$search_found_weapons_s = 
  hit_rate_traffic_aus$search_found_weapons_s / 
  max(hit_rate_traffic_aus$search_found_weapons_s, na.rm = TRUE)
hit_rate_traffic_aus$search_found_substances_s = 
  hit_rate_traffic_aus$search_found_substances_s / 
  max(hit_rate_traffic_aus$search_found_substances_s, na.rm = TRUE)

corcontaus1 = with(hit_rate_traffic_aus, lm_robust(search_found_cash_m ~ search_found_substances_m,
                                     fixed_effects = ~ monthfe))
corcontaus2 = with(hit_rate_traffic_aus, lm_robust(search_found_weapons_m ~ search_found_substances_m,
                                     fixed_effects = ~ monthfe))
corcontaus3 = with(hit_rate_traffic_aus, lm_robust(search_found_weapons_m ~ search_found_cash_m,
                                     fixed_effects = ~ monthfe))

corcontaus4 = with(hit_rate_traffic_aus, lm_robust(search_found_cash_s ~ search_found_substances_s,
                                     fixed_effects = ~ monthfe))
corcontaus5 = with(hit_rate_traffic_aus, lm_robust(search_found_weapons_s ~ search_found_substances_s,
                                     fixed_effects = ~ monthfe))
corcontaus6 = with(hit_rate_traffic_aus, lm_robust(search_found_weapons_s ~ search_found_cash_s,
                                     fixed_effects = ~ monthfe))


texreg(l = list(corcontaus1, corcontaus2, corcontaus3, corcontaus4, corcontaus5, corcontaus6),
       include.ci = FALSE,
       include.rmse = FALSE,
       include.adjrs = FALSE,
       custom.coef.map = list('search_found_substances_m' = "Drugs (Mean)",
                              'search_found_cash_m' = "Cash (Mean)",
                              'search_found_substances_s' = "Drugs (Sum)",
                              'search_found_cash_s' = "Cash (Sum)"),
       custom.model.names = c("Cash (Mean)", "Weapons (Mean)", "Weapons (Mean)",
                              "Cash (Sum)", "Weapons (Sum)", "Weapons (Sum)"),
       float.pos = "!htbp",
       booktabs = TRUE,
       caption.above = TRUE,
       caption = "\\textbf{Different types of hit rates are correlated with each other (Austin)}",
       label = "table:aushitcorr",
       scalebox = .7)


# LA stop data 

load("data/stopdata_la.RData")

# ped hit rate

hit_rate_ped_la = stopdat1 %>% 
  filter(reason_for_stop != "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            ped_contraband = mean(contraband_recovered, na.rm = TRUE),
            ped_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE)) 

hit_rate_ped_la$blmrv = 
  hit_rate_ped_la$date - min(hit_rate_ped_la$date[hit_rate_ped_la$blm == 1])

# vehicle hit rate

hit_rate_veh_la = stopdat1 %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(contraband_recovered, na.rm = TRUE),
            veh_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            veh_contraband_nalc = mean(contraband_recovered_nalc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest, na.rm = TRUE),
            arrest_felony_m = mean(arrest_felony, na.rm = TRUE),
            arrest_infraction_m = mean(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_m = mean(arrest_misdemeanor, na.rm = TRUE),
            arrest_felony_s = sum(arrest_felony, na.rm = TRUE),
            arrest_infraction_s = sum(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_s = sum(arrest_misdemeanor, na.rm = TRUE),
            contraband_recovered_evidence_s = sum(contraband_recovered_evidence, na.rm = TRUE),
            contraband_recovered_evidence_m = mean(contraband_recovered_evidence, na.rm = TRUE),
            contraband_recovered_substances_s = sum(contraband_recovered_substances, na.rm = TRUE),
            contraband_recovered_substances_m = mean(contraband_recovered_substances, na.rm = TRUE),
            contraband_recovered_weapons_s = sum(contraband_recovered_weapons, na.rm = TRUE),
            contraband_recovered_weapons_m = mean(contraband_recovered_weapons, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_veh_la$blmrv = 
  hit_rate_veh_la$date - min(hit_rate_veh_la$date[hit_rate_veh_la$blm == 1])

# assessing correlation between contraband recovery types

hit_rate_veh_la[, c("contraband_recovered_substances_m",
                    'contraband_recovered_weapons_m',
                    'contraband_recovered_evidence_m',
                    "contraband_recovered_substances_s",
                    'contraband_recovered_weapons_s',
                    'contraband_recovered_evidence_s')]  = 
  hit_rate_veh_la[, c("contraband_recovered_substances_m",
                    'contraband_recovered_weapons_m',
                    'contraband_recovered_evidence_m',
                    "contraband_recovered_substances_s",
                    'contraband_recovered_weapons_s',
                    'contraband_recovered_evidence_s')] %>% 
  apply(X = ., MARGIN = 2, FUN = function(x) x / max(x, na.rm = TRUE)) %>% 
  data.frame

hit_rate_veh_la$month = substring(hit_rate_veh_la$date, 6, 7)

corcont1 = 
  lm_robust(contraband_recovered_substances_m ~ contraband_recovered_weapons_m,
            data = hit_rate_veh_la,
            fixed_effects = ~ month)
corcont2 = 
  lm_robust(contraband_recovered_evidence_m ~ contraband_recovered_weapons_m,
            data = hit_rate_veh_la,
            fixed_effects = ~ month)
corcont3 = 
  lm_robust(contraband_recovered_evidence_m ~ contraband_recovered_substances_m,
            data = hit_rate_veh_la,
            fixed_effects = ~ month)

corcont4 = 
  lm_robust(contraband_recovered_substances_s ~ contraband_recovered_weapons_s,
            data = hit_rate_veh_la,
            fixed_effects = ~ month)
corcont5 = 
  lm_robust(contraband_recovered_evidence_s ~ contraband_recovered_weapons_s,
            data = hit_rate_veh_la,
            fixed_effects = ~ month)
corcont6 = 
  lm_robust(contraband_recovered_evidence_s ~ contraband_recovered_substances_s,
            data = hit_rate_veh_la,fixed_effects = ~ month)

texreg(l = list(corcont1, corcont2, corcont3, corcont4, corcont5, corcont6),
       include.ci = FALSE,
       include.rmse = FALSE,
       include.adjrs = FALSE,
       custom.coef.map = list('contraband_recovered_weapons_m' = "Weapons (Mean)",
                              'contraband_recovered_substances_m' = "Substances (Mean)",
                              'contraband_recovered_weapons_s' = "Weapons (Sum)",
                              'contraband_recovered_substances_s' = "Substances (Sum)"),
       custom.model.names = c("Substances (Mean)", "Evidence (Mean)", "Evidence (Mean)",
                              "Substances (Sum)", "Evidence (Sum)", "Evidence (Sum)"),
       float.pos = "!htbp",
       booktabs = TRUE,
       caption.above = TRUE,
       caption = "\\textbf{Different types of hit rates are correlated with each other (Los Angeles)}",
       label = "table:lahitcorr",
       scalebox = .7)

# consistent hit rate measures

hit_rate_veh_la_alt = stopdat1 %>% 
  filter(!(offense_code_arrest %in% c(22134, 22123, 66164, 65001, 50112,
                                      00003, 65000, 28024, 29071))) %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(contraband_recovered, na.rm = TRUE),
            veh_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            veh_contraband_nalc = mean(contraband_recovered_nalc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest, na.rm = TRUE),
            arrest_felony_m = mean(arrest_felony, na.rm = TRUE),
            arrest_infraction_m = mean(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_m = mean(arrest_misdemeanor, na.rm = TRUE),
            arrest_felony_s = sum(arrest_felony, na.rm = TRUE),
            arrest_infraction_s = sum(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_s = sum(arrest_misdemeanor, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_veh_la$blmrv = 
  hit_rate_veh_la$date - min(hit_rate_veh_la$date[hit_rate_veh_la$blm == 1])
hit_rate_veh_la_alt$blmrv = 
  hit_rate_veh_la_alt$date - min(hit_rate_veh_la_alt$date[hit_rate_veh_la_alt$blm == 1])

# Philly stops 

load(file = "data/ph_stops_ped.RData")
load(file = "data/ph_stops_veh.RData")

# hit rates (pedestrian)

hit_rate_ped_phl = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            ped_contraband = mean(ped_contraband, na.rm = TRUE),
            arrest = mean(individual_arrested, na.rm = TRUE),
            arrest_sum = sum(individual_arrested, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

hit_rate_ped_phl$blmrv = 
  hit_rate_ped_phl$date - min(hit_rate_ped_phl$date[hit_rate_ped_phl$blm == 1])

# hit rates (vehicle)

hit_rate_veh_phl = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(veh_contraband, na.rm = TRUE),
            arrest = mean(individual_arrested, na.rm = TRUE),
            arrest_sum = sum(individual_arrested, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

hit_rate_veh_phl$blmrv = 
  hit_rate_veh_phl$date - min(hit_rate_veh_phl$date[hit_rate_veh_phl$blm == 1])

# Seattle Terry Stop Data

load("data/terry.RData")

terry = terry %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE", 1, 0))

hit_rate_sea = terry %>% 
  mutate(count = 1) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  mutate(hit = ifelse(stop_resolution == "Arrest" | 
                        stop_resolution == "Citation / Infraction" | 
                        stop_resolution == "Offense Report" | 
                        stop_resolution == "Referred for Prosecution", 1, 0),
         citation = ifelse(stop_resolution == "Citation / Infraction", 1, 0),
         offense_report = ifelse(stop_resolution == "Offense Report", 1, 0),
         prosec_referral = ifelse(stop_resolution == "Referred for Prosecution", 1, 0)) %>% 
  summarize(hit = mean(hit, na.rm = TRUE),
            count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest2, na.rm = TRUE),
            arrest_sum = sum(arrest2, na.rm = TRUE),
            citation = mean(citation, na.rm = TRUE),
            offense_report = mean(offense_report, na.rm = TRUE),
            prosec_referral = mean(prosec_referral, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_sea$blmrv = 
  hit_rate_sea$date - min(hit_rate_sea$date[hit_rate_sea$blm == 1])

#### analysis: hit rates #### 

m_aus_hit_traffic = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = hit,
                                      kernel = "uni", p = 1, all = TRUE))

m_la_hit_ped = 
  with(hit_rate_ped_la, rdrobust(x = blmrv, y = ped_contraband,
                                 kernel = "uni", p = 1, all = TRUE))

m_la_hit_traffic = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband,
                                 kernel = "uni", p = 1, all = TRUE))

m_phl_hit_ped = 
  with(hit_rate_ped_phl, rdrobust(x = blmrv, y = ped_contraband,
                                  kernel = "uni", p = 1, all = TRUE))

m_phl_hit_traffic = 
  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = veh_contraband,
                                  kernel = "uni", p = 1, all = TRUE))

m_sea_hit = 
  with(hit_rate_sea, rdrobust(x = blmrv, y = hit,
                              kernel = "uni", p = 1, all = TRUE))

df_hit_rate_plot = 
  data.frame(
    
    est = c(m_aus_hit_traffic$coef[1],
            m_la_hit_ped$coef[1],
            m_la_hit_traffic$coef[1],
            m_phl_hit_ped$coef[1],
            m_phl_hit_traffic$coef[1],
            m_sea_hit$coef[1]),
    
    se = c(m_aus_hit_traffic$se[3],
           m_la_hit_ped$se[3],
           m_la_hit_traffic$se[3],
           m_phl_hit_ped$se[3],
           m_phl_hit_traffic$se[3],
           m_sea_hit$se[3]),
    
    pv = c(m_aus_hit_traffic$pv[3],
           m_la_hit_ped$pv[3],
           m_la_hit_traffic$pv[3],
           m_phl_hit_ped$pv[3],
           m_phl_hit_traffic$pv[3],
           m_sea_hit$pv[3]),
    
    city = factor(c("AUS", "LA", "LA", "PHL", "PHL", "SEA"),
                  levels = c("AUS", "LA", "PHL", "SEA")),
    
    outcome = factor(c("Hit Rate (Traffic)", "Hit Rate (Pedestrian)",
                       "Hit Rate (Traffic)", "Hit Rate (Pedestrian)",
                       "Hit Rate (Traffic)", "Hit Rate (Terry)"),
                     levels = c("Hit Rate (Traffic)", "Hit Rate (Pedestrian)",
                                "Hit Rate (Terry)"))
  )

df_hit_rate_plot = df_hit_rate_plot %>% 
  filter(!grepl(x = outcome, "estrian"))

meta_res = 
  metagen(TE = df_hit_rate_plot$est,
          seTE = df_hit_rate_plot$se, 
          sm = "MD",
          method.random.ci = 'hk')

df_hit_rate_plot = 
  bind_rows(df_hit_rate_plot, data.frame(
    est = meta_res$TE.random,
    se = meta_res$seTE.random,
    pv = meta_res$pval.random,
    city = "Meta-\nAnalysis",
    outcome = "Meta-Analysis"
  )) %>% 
  mutate(city = factor(city,
                       levels = c("AUS", "LA", "PHL", "SEA", "Meta-\nAnalysis")),
         outcome = factor(outcome,
                          levels = c("Hit Rate (Traffic)", "Hit Rate (Terry)",
                                     "Meta-Analysis")))

segment1 = df_hit_rate_plot %>% 
  ggplot() + 
  geom_hline(yintercept = 0,
             size = .2) + 
  geom_point(aes(x = city,
                 y = est,
                 shape = outcome),
             position = position_dodge(.25),
             size = 1.25) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .4) +
  labs(x = "City", y = "Coefficient\n(BLM Protest)",
       title = "A. Hit Rates",
       shape = "Outcome") + 
  theme_tufte(base_size = 9) +
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE))

#### analysis: hit rates using similar measures across austin and los angeles #### 

m_aus_hit_traffic = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = hit,
                                      kernel = "uni", p = 1, all = TRUE))

m_la_hit_traffic = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                 kernel = "uni", p = 1, all = TRUE))

m_la_hit_traffic_nalc = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_nalc,
                                 kernel = "uni", p = 1, all = TRUE))

df_hit_rate_plot = 
  data.frame(
    
    est = c(m_aus_hit_traffic$coef[1],
            m_la_hit_ped$coef[1],
            m_la_hit_traffic$coef[1],
            m_phl_hit_ped$coef[1],
            m_phl_hit_traffic$coef[1],
            m_sea_hit$coef[1]),
    
    se = c(m_aus_hit_traffic$se[3],
           m_la_hit_ped$se[3],
           m_la_hit_traffic$se[3],
           m_phl_hit_ped$se[3],
           m_phl_hit_traffic$se[3],
           m_sea_hit$se[3]),
    
    pv = c(m_aus_hit_traffic$pv[3],
           m_la_hit_ped$pv[3],
           m_la_hit_traffic$pv[3],
           m_phl_hit_ped$pv[3],
           m_phl_hit_traffic$pv[3],
           m_sea_hit$pv[3]),
    
    city = factor(c("AUS", "LA", "LA", "PHL", "PHL", "SEA"),
                  levels = c("AUS", "LA", "PHL", "SEA")),
    
    outcome = factor(c("Hit Rate (Traffic)", "Hit Rate (Pedestrian)",
                       "Hit Rate (Traffic)", "Hit Rate (Pedestrian)",
                       "Hit Rate (Traffic)", "Hit Rate (Terry)"),
                     levels = c("Hit Rate (Traffic)", "Hit Rate (Pedestrian)",
                                "Hit Rate (Terry)"))
  )

df_hit_rate_plot = df_hit_rate_plot %>% 
  filter(!grepl(x = outcome, "estrian"))

meta_res = 
  metagen(TE = df_hit_rate_plot$est,
          seTE = df_hit_rate_plot$se, 
          sm = "MD")

df_hit_rate_plot = 
  bind_rows(df_hit_rate_plot, data.frame(
    est = meta_res$TE.random,
    se = meta_res$seTE.random,
    pv = meta_res$pval.random,
    city = "Meta-\nAnalysis",
    outcome = "Meta-Analysis"
  )) %>% 
  mutate(city = factor(city,
                       levels = c("AUS", "LA", "PHL", "SEA", "Meta-\nAnalysis")),
         outcome = factor(outcome,
                          levels = c("Hit Rate (Traffic)", "Hit Rate (Terry)",
                                     "Meta-Analysis")))

segment_robust = 
  df_hit_rate_plot %>% 
  ggplot() + 
  geom_hline(yintercept = 0,
             size = .2) + 
  geom_point(aes(x = city,
                 y = est,
                 shape = outcome),
             position = position_dodge(.25),
             size = 1.25) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .4) +
  labs(x = "City", y = "Coefficient\n(BLM Protest)",
       title = "A. Hit Rates",
       shape = "Outcome") + 
  theme_tufte(base_size = 9) +
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE))

# is the LA negative effect robust? 

out1cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                 kernel = "uni", p = 2, all = TRUE))
out2cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "uni", p = 3, all = TRUE))
out3cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "tri", p = 1, all = TRUE))
out4cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "tri", p = 2, all = TRUE))
out5cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "tri", p = 3, all = TRUE))
out6cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "epa", p = 1, all = TRUE))
out7cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "epa", p = 2, all = TRUE))
out8cla = with(hit_rate_veh_la, rdrobust(x = blmrv, y = veh_contraband_alc,
                                         kernel = "epa", p = 3, all = TRUE))

butsegment_notrobust = data.frame(
  
  est = c(out1cla$coef[3],
    out2cla$coef[3],
    out3cla$coef[3],
    out4cla$coef[3],
    out5cla$coef[3],
    out6cla$coef[3],
    out7cla$coef[3],
    out8cla$coef[3]),
  
  se = c(out1cla$se[3],
    out2cla$se[3],
    out3cla$se[3],
    out4cla$se[3],
    out5cla$se[3],
    out6cla$se[3],
    out7cla$se[3],
    out8cla$se[3]),
  
  pv = c(out1cla$pv[3],
         out2cla$pv[3],
         out3cla$pv[3],
         out4cla$pv[3],
         out5cla$pv[3],
         out6cla$pv[3],
         out7cla$pv[3],
         out8cla$pv[3]),
  
  poly = c(2, 3, rep(c(1, 2, 3), 2)),
  
  kernel = c('Uni.', "Uni.", rep('Tri.', 3),  rep('Epa.', 3))
  
  
) %>% 
  mutate(poly_kern = paste0(kernel, ',\n', poly)) %>% 
  mutate(poly_kern = factor(poly_kern, levels = c("Uni.,\n2", "Uni.,\n3",
                                                  "Tri.,\n1", "Tri.,\n2",
                                                  "Tri.,\n3", "Epa.,\n1",
                                                  "Epa.,\n2", "Epa.,\n3"))) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = poly_kern, y = est),
             size = 2) + 
  geom_errorbar(aes(x = poly_kern, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                width = 0,
                size = .8) + 
  geom_errorbar(aes(x = poly_kern, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                width = 0,
                size = .6) + 
  labs(x = "Polynomial, Kernel", y = "RDiT Coefficient (BLM Protest)",
       title = "B. Assessing Robustness of LA Effect") + 
  theme_tufte(base_size = 9)

segment_robust_grob =
  arrangeGrob(segment_robust, butsegment_notrobust, ncol = 2)

ggsave(plot = segment_robust_grob,
       filename = "pics/robconsistentmeasure.jpeg",
       dpi = 300, width = 6, height = 2.5)

#### analysis: arrest rates #### 

m_aus_arr_traffic = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = arrest,
                                      kernel = "uni", p = 1, all = TRUE))

m_la_arr_ped = 
  with(hit_rate_ped_la, rdrobust(x = blmrv, y = arrest,
                                 kernel = "uni", p = 1, all = TRUE))

m_la_arr_traffic = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = arrest,
                                 kernel = "uni", p = 1, all = TRUE))

m_phl_arr_ped = 
  with(hit_rate_ped_phl, rdrobust(x = blmrv, y = arrest,
                                  kernel = "uni", p = 1, all = TRUE))

m_phl_arr_traffic = 
  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = arrest,
                                  kernel = "uni", p = 1, all = TRUE))
m_sea_arr = 
  with(subset(hit_rate_sea, date >= as.Date("2019-06-01")),
       rdrobust(x = blmrv, y = arrest,
                              kernel = "uni", p = 1, all = TRUE))

## assessing sd outcome estimates
#
#m_aus_arr_traffic = 
#  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = scale(arrest),
#                                      kernel = "uni", p = 1, all = TRUE))
#
#m_la_arr_ped = 
#  with(hit_rate_ped_la, rdrobust(x = blmrv, y = scale(arrest),
#                                 kernel = "uni", p = 1, all = TRUE))
#
#m_la_arr_traffic = 
#  with(hit_rate_veh_la, rdrobust(x = blmrv, y = scale(arrest),
#                                 kernel = "uni", p = 1, all = TRUE))
#
#m_phl_arr_ped = 
#  with(hit_rate_ped_phl, rdrobust(x = blmrv, y = scale(arrest),
#                                  kernel = "uni", p = 1, all = TRUE))
#
#m_phl_arr_traffic = 
#  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = scale(arrest),
#                                  kernel = "uni", p = 1, all = TRUE))
#m_sea_arr = 
#  with(subset(hit_rate_sea, date >= as.Date("2019-06-01")),
#       rdrobust(x = blmrv, y = scale(arrest),
#                kernel = "uni", p = 1, all = TRUE))
#
#m_aus_arr_traffic$coef[1] * 100
#m_la_arr_ped$coef[1] * 100
#m_la_arr_traffic$coef[1] * 100
#m_phl_arr_ped$coef[1] * 100
#m_phl_arr_traffic$coef[1] * 100
#m_sea_arr$coef[1] * 100


# stuff

df_arr_rate_plot = 
  data.frame(
    
    est = c(m_aus_arr_traffic$coef[1],
            m_la_arr_ped$coef[1],
            m_la_arr_traffic$coef[1],
            m_phl_arr_ped$coef[1],
            m_phl_arr_traffic$coef[1],
            m_sea_arr$coef[1]),
    
    se = c(m_aus_arr_traffic$se[3],
           m_la_arr_ped$se[3],
           m_la_arr_traffic$se[3],
           m_phl_arr_ped$se[3],
           m_phl_arr_traffic$se[3],
           m_sea_arr$se[3]),
    
    pv = c(m_aus_arr_traffic$pv[3],
           m_la_arr_ped$pv[3],
           m_la_arr_traffic$pv[3],
           m_phl_arr_ped$pv[3],
           m_phl_arr_traffic$pv[3],
           m_sea_arr$pv[3]),
    
    city = factor(c("AUS", "LA", "LA", "PHL", "PHL", "SEA"),
                  levels = c("AUS", "LA", "PHL", "SEA")),
    
    outcome = factor(c("Arrest Rate (Traffic)", "Arrest Rate (Pedestrian)",
                       "Arrest Rate (Traffic)", "Arrest Rate (Pedestrian)",
                       "Arrest Rate (Traffic)", "Arrest Rate (Terry)"),
                     levels = c("Arrest Rate (Traffic)", "Arrest Rate (Pedestrian)",
                                "Arrest Rate (Terry)"))
  )


df_arr_rate_plot = 
  df_arr_rate_plot %>% filter(!grepl(x = outcome, "estrian")) 

meta_res = 
  metagen(TE = df_arr_rate_plot$est,
          seTE = df_arr_rate_plot$se, 
          sm = "MD",
          method.random.ci = "HK")

df_arr_rate_plot = 
  bind_rows(df_arr_rate_plot, data.frame(
    est = meta_res$TE.random,
    se = meta_res$seTE.random,
    pv = meta_res$pval.random,
    city = "Meta-\nAnalysis",
    outcome = "Meta-Analysis"
  )) %>% 
  mutate(city = factor(city,
                       levels = c("AUS", "LA", "PHL",
                                  "SEA", "Meta-\nAnalysis")),
         outcome = factor(outcome,
                          levels = c("Arrest Rate (Traffic)",
                                     "Arrest Rate (Terry)",
                                     "Meta-Analysis")))

segment2 = df_arr_rate_plot %>% 
  ggplot() + 
  geom_hline(yintercept = 0,
             size = .2) + 
  geom_point(aes(x = city,
                 y = est,
                 shape = outcome),
             position = position_dodge(.25),
             size = 1.25) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .4) +
  labs(x = "City", y = "Coefficient\n(BLM Protest)",
       title = "B. Arrest Rates",
       shape = "Outcome") + 
  theme_tufte(base_size = 9) +
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE))

#### arrest rates, placebos before and after #### 

list_to_loop_df = 
  list(hit_rate_traffic_aus, hit_rate_veh_la, hit_rate_veh_phl, hit_rate_sea)
city_label = c("AUS", "LA", "PHL", "SEA")
list_to_loop_df_out = as.list(rep(NA, length(list_to_loop_df)))
seq_2_l = seq(from = -100, to = 100, by = 20)
stop_label = c("Traffic", "Vehicle", "Vehicle", "Terry")

for (i in 1:length(list_to_loop_df)) {
  
  print(paste0("I iteration ", i))
  
  df_seq_out = 
    matrix(NA, nrow = length(seq_2_l), ncol = 6) %>% 
    as.data.frame %>% 
    `colnames<-` (c("est", "se", "pv", 'plcval', 'city', 'outcome'))
  
  for (k in 1:length(seq_2_l)) {
    
    print(paste0("K iteration ", k))
    
    out = with(list_to_loop_df[[i]],
               rdrobust(x = I(blmrv + seq_2_l[k]), y = arrest, p = 1,
                        kernel = 'uni'))
    
    df_seq_out[k, 1] = out$coef[3]
    df_seq_out[k, 2] = out$se[3]
    df_seq_out[k, 3] = out$pv[3]
    df_seq_out[k, 4] = seq_2_l[k]
    df_seq_out[k, 5] = city_label[i]
    df_seq_out[k, 6] = stop_label[i]
    
  }
  
  list_to_loop_df_out[[i]] = df_seq_out
  
}

list_to_loop_df_out %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(city_stop = paste0(city, ', ', outcome)) %>% 
  mutate(plcvalbin = ifelse(plcval == 0, 1, 0)) %>% 
  mutate(plcvalbin = factor(plcvalbin, levels = c(0, 1))) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = plcval, y = est,
                 col = plcvalbin),
             size = 2) + 
  geom_errorbar(aes(x = plcval, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = plcvalbin),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = plcval, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = plcvalbin),
                size = .8,
                width = 0) + 
  guides(col = 'none') + 
  labs(x = "Placebo Discontinuity (Days to BLM Protest)",
       y = "RDiT Coefficient",
       title = 'Placebo Discontinuities Before and After the BLM Protest (Arrest Rate Outcome)') + 
  scale_color_grey(start = .6, end = 0) + 
  facet_wrap(~city_stop, scales = 'free') + 
  theme_tufte()

ggsave(plot = last_plot(), filename = 'pics/plcdisc_arr.jpeg', dpi = 600,
       width = 8, height = 4.5)

#### analysis: mechanism test on arrests #### 

# effect on counts of stops
# effect on counts of arrests
# stop count effect relative to baseline average stop number
# arrest count effect relative to baseline average arrest number

# austin reconversion 

arr_avg_preblm_aus = 
  mean(hit_rate_traffic_aus$arrest_sum[hit_rate_traffic_aus$date < hit_rate_traffic_aus$date[hit_rate_traffic_aus$blmrv == 0] & 
  hit_rate_traffic_aus$date >= as.Date("2020-01-01")])
cnt_avg_preblm_aus = 
  mean(hit_rate_traffic_aus$count[hit_rate_traffic_aus$date < hit_rate_traffic_aus$date[hit_rate_traffic_aus$blmrv == 0] & 
                                         hit_rate_traffic_aus$date >= as.Date("2020-01-01")])
hit_rate_traffic_aus$count_sum_ratio =
  hit_rate_traffic_aus$count / cnt_avg_preblm_aus
hit_rate_traffic_aus$arrest_sum_ratio = 
  hit_rate_traffic_aus$arrest_sum / arr_avg_preblm_aus 

# los angeles reconversion 

arr_avg_preblm_la = 
  mean(hit_rate_veh_la$arrest_sum[hit_rate_veh_la$date < hit_rate_veh_la$date[hit_rate_veh_la$blmrv == 0] & 
                                    hit_rate_veh_la$date >= as.Date("2020-01-01")])
cnt_avg_preblm_la = 
  mean(hit_rate_traffic_aus$count[hit_rate_veh_la$date < hit_rate_veh_la$date[hit_rate_veh_la$blmrv == 0] & 
                                    hit_rate_veh_la$date >= as.Date("2020-01-01")])

hit_rate_veh_la$count_sum_ratio =
  hit_rate_veh_la$count / cnt_avg_preblm_la
hit_rate_veh_la$arrest_sum_ratio = 
  hit_rate_veh_la$arrest_sum / arr_avg_preblm_la

# philadelphia reconversion 

arr_avg_preblm_ph = 
  mean(hit_rate_veh_phl$arrest_sum[hit_rate_veh_phl$date < hit_rate_veh_phl$date[hit_rate_veh_phl$blmrv == 0] & 
                                     hit_rate_veh_phl$date >= as.Date("2020-01-01")])
cnt_avg_preblm_ph = 
  mean(hit_rate_veh_phl$count[hit_rate_veh_phl$date < hit_rate_veh_phl$date[hit_rate_veh_phl$blmrv == 0] & 
                                hit_rate_veh_phl$date >= as.Date("2020-01-01")])

hit_rate_veh_phl$count_sum_ratio =
  hit_rate_veh_phl$count / cnt_avg_preblm_ph
hit_rate_veh_phl$arrest_sum_ratio = 
  hit_rate_veh_phl$arrest_sum / arr_avg_preblm_ph

# seattle reconversion 

arr_avg_preblm_sea = 
  mean(hit_rate_sea$arrest_sum[hit_rate_sea$date < hit_rate_sea$date[hit_rate_sea$blmrv == 0] & 
                                     hit_rate_sea$date >= as.Date("2020-01-01")])
cnt_avg_preblm_sea = 
  mean(hit_rate_sea$count[hit_rate_sea$date < hit_rate_sea$date[hit_rate_sea$blmrv == 0] & 
                            hit_rate_sea$date >= as.Date("2020-01-01")])

hit_rate_sea$count_sum_ratio =
  hit_rate_sea$count / cnt_avg_preblm_sea
hit_rate_sea$arrest_sum_ratio = 
  hit_rate_sea$arrest_sum / arr_avg_preblm_sea

# estimates - austin

aus_arr_mech_num = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = arrest_sum,
                                      kernel = "uni", p = 1, all = TRUE))
aus_arr_mech_denom = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = count,
                                      kernel = "uni", p = 1, all = TRUE))
aus_arr_mech_num_ratio = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = arrest_sum_ratio,
                                      kernel = "uni", p = 1, all = TRUE))
aus_arr_mech_denom_ratio = 
  with(hit_rate_traffic_aus, rdrobust(x = blmrv, y = count_sum_ratio,
                                      kernel = "uni", p = 1, all = TRUE))

# estimates - los angeles

la_arr_mech_num = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = arrest_sum,
                                      kernel = "uni", p = 1, all = TRUE))
la_arr_mech_denom = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = count,
                                      kernel = "uni", p = 1, all = TRUE))
la_arr_mech_num_ratio = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = arrest_sum_ratio,
                                      kernel = "uni", p = 1, all = TRUE))
la_arr_mech_denom_ratio = 
  with(hit_rate_veh_la, rdrobust(x = blmrv, y = count_sum_ratio,
                                      kernel = "uni", p = 1, all = TRUE))

hit_rate_veh_la_alt$arrest_felony_m2 = 
  hit_rate_veh_la_alt$arrest_felony_s / hit_rate_veh_la_alt$count
hit_rate_veh_la_alt$arrest_misdemeanor_m2 = 
  hit_rate_veh_la_alt$arrest_misdemeanor_s / hit_rate_veh_la_alt$count
hit_rate_veh_la_alt$arrest_infraction_s = 
  hit_rate_veh_la_alt$arrest_infraction_s / hit_rate_veh_la_alt$count


# estimates - philadelphia

phl_arr_mech_num = 
  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = arrest_sum,
                               kernel = "uni", p = 1, all = TRUE))
phl_arr_mech_denom = 
  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = count,
                                 kernel = "uni", p = 1, all = TRUE))
phl_arr_mech_num_ratio = 
  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = arrest_sum_ratio,
                                 kernel = "uni", p = 1, all = TRUE))
phl_arr_mech_denom_ratio = 
  with(hit_rate_veh_phl, rdrobust(x = blmrv, y = count_sum_ratio,
                                 kernel = "uni", p = 1, all = TRUE))

# estimates - seattle

sea_arr_mech_num = 
  with(subset(hit_rate_sea, date >= as.Date("2019-06-01")),
       rdrobust(x = blmrv, y = arrest_sum,
                                  kernel = "uni", p = 1, all = TRUE))
sea_arr_mech_denom = 
  with(subset(hit_rate_sea, date >= as.Date("2019-06-01")),
       rdrobust(x = blmrv, y = count,
                                  kernel = "uni", p = 1, all = TRUE))
sea_arr_mech_num_ratio = 
  with(subset(hit_rate_sea, date >= as.Date("2019-06-01")),
       rdrobust(x = blmrv, y = arrest_sum_ratio,
                                  kernel = "uni", p = 1, all = TRUE))
sea_arr_mech_denom_ratio = 
  with(subset(hit_rate_sea, date >= as.Date("2019-06-01")),
       rdrobust(x = blmrv, y = count_sum_ratio,
                                  kernel = "uni", p = 1, all = TRUE))

# now, plotting....
# panel A: arrest counts by city
# panel B: stop counts by city
# panel C: stop vs arrest ratio (Austin)
# panel D: stop vs arrest ratio (LA)
# panel E: stop vs arrest ratio (philly)
# panel F: stop vs arrest ratio (seattle)

discretp1 = data.frame(
  
  est = 
    c(aus_arr_mech_num$coef[3], 
      la_arr_mech_num$coef[3], 
      phl_arr_mech_num$coef[3], 
      sea_arr_mech_num$coef[3]),
  
  se = 
    c(aus_arr_mech_num$se[3],
    la_arr_mech_num$se[3],
    phl_arr_mech_num$se[3], 
    sea_arr_mech_num$se[3]),
  
  pv = 
    c(aus_arr_mech_num$pv[3],
      la_arr_mech_num$pv[3],
      phl_arr_mech_num$pv[3], 
      sea_arr_mech_num$pv[3]),
  
  city = c("AUS", "LA", "PHL", "SEA"),
  
  outcome = c("# Arrests (Traffic)",
              "# Arrests (Traffic)",
              "# Arrests (Traffic)",
              "# Arrests (Terry)")
  
) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est, shape = outcome),
             size = 2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                size = .8,
                width = 0) + 
  labs(x = "City",
       y = "RDiT Coefficient\n(BLM Protest)",
       title = "A. Total Arrests (Numerator)",
       shape = "Outcome") + 
  theme_tufte() + 
  theme(legend.key.size = unit(.1, 'cm'),
      legend.title = element_text(size = 6),
      legend.text = element_text(size = 4),
      legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 1, byrow = TRUE))

discretp2 = data.frame(
  
  est = 
    c(aus_arr_mech_denom$coef[3], 
      la_arr_mech_denom$coef[3], 
      phl_arr_mech_denom$coef[3], 
      sea_arr_mech_denom$coef[3]),
  
  se = 
    c(aus_arr_mech_denom$se[3],
      la_arr_mech_denom$se[3],
      phl_arr_mech_denom$se[3], 
      sea_arr_mech_denom$se[3]),
  
  pv = 
    c(aus_arr_mech_denom$pv[3],
      la_arr_mech_denom$pv[3],
      phl_arr_mech_denom$pv[3], 
      sea_arr_mech_denom$pv[3]),
  
  city = c("AUS", "LA", "PHL", "SEA"),
  
  outcome = c("# Stops (Traffic)",
              "# Stops (Traffic)",
              "# Stops (Traffic)",
              "# Stops (Terry)")
  
) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est, shape = outcome),
             size = 2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                size = .8,
                width = 0) + 
  labs(x = "City",
       y = "RDiT Coefficient\n(BLM Protest)",
       title = "B. Total Stops (Denominator)",
       shape = "Outcome") + 
  theme_tufte() + 
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 1, byrow = TRUE))

ratiosresdf = data.frame(
  
  est = 
    c(aus_arr_mech_num_ratio$coef[3], 
      la_arr_mech_num_ratio$coef[3], 
      phl_arr_mech_num_ratio$coef[3], 
      sea_arr_mech_num_ratio$coef[3],
      aus_arr_mech_denom_ratio$coef[3], 
      la_arr_mech_denom_ratio$coef[3], 
      phl_arr_mech_denom_ratio$coef[3], 
      sea_arr_mech_denom_ratio$coef[3]),
  
  se = 
    c(aus_arr_mech_num_ratio$se[3],
      la_arr_mech_num_ratio$se[3],
      phl_arr_mech_num_ratio$se[3], 
      sea_arr_mech_num_ratio$se[3],
      aus_arr_mech_denom_ratio$se[3],
      la_arr_mech_denom_ratio$se[3],
      phl_arr_mech_denom_ratio$se[3], 
      sea_arr_mech_denom_ratio$se[3]),
  
  pv = 
    c(aus_arr_mech_num_ratio$pv[3],
      la_arr_mech_num_ratio$pv[3],
      phl_arr_mech_num_ratio$pv[3], 
      sea_arr_mech_num_ratio$pv[3],
      aus_arr_mech_denom_ratio$pv[3],
      la_arr_mech_denom_ratio$pv[3],
      phl_arr_mech_denom_ratio$pv[3], 
      sea_arr_mech_denom_ratio$pv[3]),
  
  city = rep(c("AUS", "LA", "PHL", "SEA"), 2),
  
  outcome = c(rep("Arrest Ratio", 4),
              "Stop Ratio",
              "Stop Ratio",
              "Stop Ratio",
              "Stop Ratio")
  
)


ratiosresdf_aus = ratiosresdf %>% filter(city == "AUS")
abs((ratiosresdf_aus$est[2] - ratiosresdf_aus$est[1]) / 
  sqrt(ratiosresdf_aus$se[2]^2 + ratiosresdf_aus$se[1]^2))

# p < 0.001

discretp3 = 
  ratiosresdf %>%
  filter(city == "AUS") %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est),
             size = 2) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                size = .8,
                width = 0) + 
  annotate("text",
           label = paste0("Stop - Arrest Ratio\nCoef. Difference: ",
                          round((ratiosresdf_aus$est[2] - 
                                   ratiosresdf_aus$est[1]), 2),
                          "\np < 0.001"),
           x = 1.5,
           y = -1,
           family = 'serif',
           size = 3
  ) + 
  labs(x = "Outcome",
       y = "RDiT Coefficient\n(BLM Protest)",
       title = "C. Arrest vs Stop Ratio (Austin)") + 
  theme_tufte() + 
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 1, byrow = TRUE))

ratiosresdf_la = ratiosresdf %>% filter(city == "LA")
abs((ratiosresdf_la$est[2] - ratiosresdf_la$est[1]) / 
      sqrt(ratiosresdf_la$se[2]^2 + ratiosresdf_la$se[1]^2))

discretp4 = ratiosresdf %>%
  filter(city == "LA") %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est),
             size = 2) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                size = .8,
                width = 0) + 
  annotate("text",
           label = paste0("Stop - Arrest Ratio\nCoef. Difference: ",
                          round((ratiosresdf_la$est[2] - 
                                   ratiosresdf_la$est[1]), 2),
                          "\np < 0.001"),
           x = 1.5,
           y = -5,
           family = 'serif',
           size = 3
  ) + 
  labs(x = "Outcome",
       y = "RDiT Coefficient\n(BLM Protest)",
       title = "D. Arrest vs Stop Ratio (LA)") + 
  theme_tufte() + 
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 1, byrow = TRUE))

ratiosresdf_phl = ratiosresdf %>% filter(city == "PHL")

abs((ratiosresdf_phl$est[2] - ratiosresdf_phl$est[1]) / 
      sqrt(ratiosresdf_phl$se[2]^2 + ratiosresdf_phl$se[1]^2))

discretp5 = ratiosresdf %>%
  filter(city == "PHL") %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est),
             size = 2) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                size = .8,
                width = 0) + 
  annotate("text",
           label = paste0("Stop - Arrest Ratio\nCoef. Difference: ",
                          round((ratiosresdf_phl$est[2] - 
                                   ratiosresdf_phl$est[1]), 2),
                          "\np < 0.05"),
           x = 1.5,
           y = -.25,
           family = 'serif',
           size = 3
  ) + 
  labs(x = "Outcome",
       y = "RDiT Coefficient\n(BLM Protest)",
       title = "E. Arrest vs Stop Ratio (PHL)") + 
  theme_tufte() + 
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 1, byrow = TRUE))

ratiosresdf_sea = ratiosresdf %>% filter(city == "SEA")

abs((ratiosresdf_sea$est[2] - ratiosresdf_sea$est[1]) / 
      sqrt(ratiosresdf_sea$se[2]^2 + ratiosresdf_sea$se[1]^2))

discretp6 = ratiosresdf %>%
  filter(city == "SEA") %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est),
             size = 2) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = outcome, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                size = .8,
                width = 0) + 
  annotate("text",
           label = paste0("Stop - Arrest Ratio\nCoef. Difference: ",
                          round((ratiosresdf_sea$est[2] - 
                                   ratiosresdf_sea$est[1]), 2),
                          "\np = 0.06"),
           x = 1.5,
           y = -.45,
           family = 'serif',
           size = 3
  ) + 
  labs(x = "Outcome",
       y = "RDiT Coefficient\n(BLM Protest)",
       title = "F. Arrest vs Stop Ratio (SEA)") + 
  theme_tufte() + 
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 1, byrow = TRUE))

discret_grob = 
  arrangeGrob(discretp1, discretp2, discretp3, discretp4, discretp5, discretp6,
            ncol = 2)

ggsave(plot = discret_grob, filename = "pics/discret_rc.jpeg", dpi = 600, width = 8,
       height = 6.5)

#### analysis: mechanism test on arrests, decomposition #### 

# ID

(head(table(stopdat1$offense_code_arrest[stopdat1$date <= (as.Date("2020-05-28") + 15) & stopdat1$date >= (as.Date("2020-05-28"))])[
  order(table(stopdat1$offense_code_arrest[stopdat1$date <= (as.Date("2020-05-28") + 15) & stopdat1$date >= (as.Date("2020-05-28"))]),
        decreasing = TRUE)
], 30))[!(names(head(table(stopdat1$offense_code_arrest[stopdat1$date <= (as.Date("2020-05-28") + 15) & stopdat1$date >= (as.Date("2020-05-28"))])[
  order(table(stopdat1$offense_code_arrest[stopdat1$date <= (as.Date("2020-05-28") + 15) & stopdat1$date >= (as.Date("2020-05-28"))]),
        decreasing = TRUE)
], 30))) %in% 
  
  (names(head(table(stopdat1$offense_code_arrest[stopdat1$date >= (as.Date("2020-05-28") - 15) & stopdat1$date < (as.Date("2020-05-28"))])[
    order(table(stopdat1$offense_code_arrest[stopdat1$date >= (as.Date("2020-05-28") - 15) & stopdat1$date < (as.Date("2020-05-28"))]),
          decreasing = TRUE)
  ], 30)))]

hit_rate_veh_la_alt_rmp = 
  stopdat1 %>% 
  filter(!(offense_code_arrest %in% c(22134, 22123, 66164, 65001, 50112, 00003, 65000, 28024, 29071))) %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(contraband_recovered, na.rm = TRUE),
            veh_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            veh_contraband_nalc = mean(contraband_recovered_nalc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest, na.rm = TRUE),
            arrest_felony_m = mean(arrest_felony, na.rm = TRUE),
            arrest_infraction_m = mean(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_m = mean(arrest_misdemeanor, na.rm = TRUE),
            arrest_felony_s = sum(arrest_felony, na.rm = TRUE),
            arrest_infraction_s = sum(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_s = sum(arrest_misdemeanor, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_veh_la_alt = 
  stopdat1 %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(contraband_recovered, na.rm = TRUE),
            veh_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            veh_contraband_nalc = mean(contraband_recovered_nalc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest, na.rm = TRUE),
            arrest_felony_m = mean(arrest_felony, na.rm = TRUE),
            arrest_infraction_m = mean(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_m = mean(arrest_misdemeanor, na.rm = TRUE),
            arrest_felony_s = sum(arrest_felony, na.rm = TRUE),
            arrest_infraction_s = sum(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_s = sum(arrest_misdemeanor, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_veh_la_alt$blmrv = 
  hit_rate_veh_la_alt$date - 
  min(hit_rate_veh_la_alt$date[hit_rate_veh_la_alt$blm == 1])
hit_rate_veh_la_alt_rmp$blmrv = 
  hit_rate_veh_la_alt_rmp$date -
  min(hit_rate_veh_la_alt_rmp$date[hit_rate_veh_la_alt_rmp$blm == 1])

# infraction sum doesn't change, misdemeanor arrests decline, felony arrests also 
# decline 

smod1 = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_infraction_s),
       p = 1, kernel = "uni")
smod2 = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_misdemeanor_s),
       p = 1, kernel = "uni")
smod3 = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_felony_s),
       p = 1, kernel = "uni")

smod1_adj = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_infraction_s),
       p = 1, kernel = "uni")
smod2_adj = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_misdemeanor_s),
       p = 1, kernel = "uni")
smod3_adj = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_felony_s),
       p = 1, kernel = "uni")

smod1_m = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_infraction_m),
       p = 1, kernel = "uni")
smod2_m = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_misdemeanor_m),
       p = 1, kernel = "uni")
smod3_m = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_felony_m),
       p = 1, kernel = "uni")

smod1_adj_m = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_infraction_m),
       p = 1, kernel = "uni")
smod2_adj_m = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_misdemeanor_m),
       p = 1, kernel = "uni")
smod3_adj_m = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_felony_m),
       p = 1, kernel = "uni")

# estimates once more, but without normalizing only among arrests, but rather
# stops 

hit_rate_veh_la_alt_rmp = 
  stopdat1 %>% 
  filter(!(offense_code_arrest %in% 
             c(22134, 22123, 66164, 65001, 50112,
               00003, 65000, 28024, 29071))) %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  mutate(arrest_felony = ifelse(is.na(arrest_felony), 0, arrest_felony),
         arrest_infraction = ifelse(is.na(arrest_infraction), 0, arrest_infraction),
         arrest_misdemeanor = ifelse(is.na(arrest_misdemeanor), 0, arrest_misdemeanor)) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(contraband_recovered, na.rm = TRUE),
            veh_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            veh_contraband_nalc = mean(contraband_recovered_nalc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest, na.rm = TRUE),
            arrest_felony_m = mean(arrest_felony, na.rm = TRUE),
            arrest_infraction_m = mean(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_m = mean(arrest_misdemeanor, na.rm = TRUE),
            arrest_felony_s = sum(arrest_felony, na.rm = TRUE),
            arrest_infraction_s = sum(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_s = sum(arrest_misdemeanor, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_veh_la_alt = 
  stopdat1 %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  mutate(arrest_felony = ifelse(is.na(arrest_felony), 0, arrest_felony),
         arrest_infraction = ifelse(is.na(arrest_infraction), 0, arrest_infraction),
         arrest_misdemeanor = ifelse(is.na(arrest_misdemeanor), 0, arrest_misdemeanor)) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(contraband_recovered, na.rm = TRUE),
            veh_contraband_alc = mean(contraband_recovered_alc, na.rm = TRUE),
            veh_contraband_nalc = mean(contraband_recovered_nalc, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest_m = mean(arrest, na.rm = TRUE),
            arrest_sum = sum(arrest, na.rm = TRUE),
            arrest_felony_m = mean(arrest_felony, na.rm = TRUE),
            arrest_infraction_m = mean(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_m = mean(arrest_misdemeanor, na.rm = TRUE),
            arrest_felony_s = sum(arrest_felony, na.rm = TRUE),
            arrest_infraction_s = sum(arrest_infraction, na.rm = TRUE),
            arrest_misdemeanor_s = sum(arrest_misdemeanor, na.rm = TRUE)) %>% 
  rename(arrest = arrest_m)

hit_rate_veh_la_alt$blmrv = 
  hit_rate_veh_la_alt$date - 
  min(hit_rate_veh_la_alt$date[hit_rate_veh_la_alt$blm == 1])
hit_rate_veh_la_alt_rmp$blmrv = 
  hit_rate_veh_la_alt_rmp$date -
  min(hit_rate_veh_la_alt_rmp$date[hit_rate_veh_la_alt_rmp$blm == 1])


smod1b_m = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_infraction_m),
       p = 1, kernel = "uni")
smod2b_m = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_misdemeanor_m),
       p = 1, kernel = "uni")
smod3b_m = 
  with(hit_rate_veh_la_alt, rdrobust(x = blmrv, y = arrest_felony_m),
       p = 1, kernel = "uni")

smod1_adjb_m = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_infraction_m),
       p = 1, kernel = "uni")
smod2_adjb_m = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_misdemeanor_m),
       p = 1, kernel = "uni")
smod3_adjb_m = 
  with(hit_rate_veh_la_alt_rmp, rdrobust(x = blmrv, y = arrest_felony_m),
       p = 1, kernel = "uni")

# big plot

decomp_plot_sum = 
  data.frame(
  est = c(smod1$coef[1], smod2$coef[1], smod3$coef[1],
          smod1_adj$coef[1], smod2_adj$coef[1], smod3_adj$coef[1]),
  se = c(smod1$se[3], smod2$se[3], smod3$se[3],
          smod1_adj$se[3], smod2_adj$se[3], smod3_adj$se[3]),
  outcome = factor(c("# Infraction", "# Misdemeanor", "# Felony"),
                   levels =c("# Infraction", "# Misdemeanor", "# Felony")),
  dataset = factor(c("Unadjusted", "Unadjusted", "Unadjusted",
              "Protest-Adjusted", "Protest-Adjusted", "Protest-Adjusted"),
              levels = c("Unadjusted", "Protest-Adjusted"))
)

decomp_plot_mean = 
  data.frame(
    est = c(smod1_m$coef[1], smod2_m$coef[1], smod3_m$coef[1],
            smod1_adj_m$coef[1], smod2_adj_m$coef[1], smod3_adj_m$coef[1]),
    se = c(smod1_m$se[3], smod2_m$se[3], smod3_m$se[3],
           smod1_adj_m$se[3], smod2_adj_m$se[3], smod3_adj_m$se[3]),
    outcome = factor(c("Pr(Infraction)", "Pr(Misdemeanor)", "Pr(Felony)"),
                     levels =c("Pr(Infraction)", "Pr(Misdemeanor)", "Pr(Felony)")),
    dataset = factor(c("Unadjusted", "Unadjusted", "Unadjusted",
                       "Protest-Adjusted", "Protest-Adjusted", "Protest-Adjusted"),
                     levels = c("Unadjusted", "Protest-Adjusted"))
  )

decomp_plot_mean2 = 
  data.frame(
    est = c(smod1b_m$coef[1], smod2b_m$coef[1], smod3b_m$coef[1],
            smod1_adjb_m$coef[1], smod2_adjb_m$coef[1], smod3_adjb_m$coef[1]),
    se = c(smod1b_m$se[3], smod2b_m$se[3], smod3b_m$se[3],
           smod1_adjb_m$se[3], smod2_adjb_m$se[3], smod3_adjb_m$se[3]),
    outcome = factor(c("Pr(Infraction)", "Pr(Misdemeanor)", "Pr(Felony)"),
                     levels =c("Pr(Infraction)", "Pr(Misdemeanor)", "Pr(Felony)")),
    dataset = factor(c("Unadjusted", "Unadjusted", "Unadjusted",
                       "Protest-Adjusted", "Protest-Adjusted", "Protest-Adjusted"),
                     levels = c("Unadjusted", "Protest-Adjusted"))
  )


decomp_plot1 = decomp_plot_sum %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est, col = dataset),
             position = position_dodge(.5),
             size = 2) +
  geom_errorbar(aes(x = outcome,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = dataset),
                position = position_dodge(.5),
                width = 0,
                size = .6) + 
  geom_errorbar(aes(x = outcome,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = dataset),
                position = position_dodge(.5),
                width = 0,
                size = .8) + 
  labs(x = "Outcome", y = "RDiT Coefficient (BLM Protest)",
       col = "Dataset",
       title = "A. Outcome: # Arrests") + 
  scale_color_grey(start = 0, end = .4) + 
  theme_tufte(base_size = 9) + 
  theme(legend.position = c(.25, .25),
        legend.key.size = unit(.2, 'cm'),
        legend.title.align = 0.5)


decomp_plot2 = decomp_plot_mean %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est, col = dataset),
             position = position_dodge(.5),
             size = 2) +
  geom_errorbar(aes(x = outcome,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = dataset),
                position = position_dodge(.5),
                width = 0,
                size = .6) + 
  geom_errorbar(aes(x = outcome,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = dataset),
                position = position_dodge(.5),
                width = 0,
                size = .8) + 
  labs(x = "Outcome", y = "",
       col = "Dataset",
       title = "B. Outcome: Arrest Type / Arrest") + 
  guides(col = 'none') +
  scale_color_grey(start = 0, end = .4) + 
  theme_tufte(base_size = 9) + 
  theme(legend.position = c(.25, .25))

decomp_plot3 = decomp_plot_mean2 %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = outcome, y = est, col = dataset),
             position = position_dodge(.5),
             size = 2) +
  geom_errorbar(aes(x = outcome,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = dataset),
                position = position_dodge(.5),
                width = 0,
                size = .6) + 
  geom_errorbar(aes(x = outcome,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = dataset),
                position = position_dodge(.5),
                width = 0,
                size = .8) + 
  labs(x = "Outcome", y = "",
       col = "Dataset",
       title = "C. Outcome: Arrest Type / Stop") + 
  guides(col = 'none') + 
  scale_color_grey(start = 0, end = .4) + 
  theme_tufte(base_size = 9) + 
  theme(legend.position = c(.25, .45))

decomp_plot_grob = 
  arrangeGrob(decomp_plot1, decomp_plot2, decomp_plot3, ncol = 3)

ggsave(plot = decomp_plot_grob, 
       filename = "pics/decomp_plot.jpeg",
       dpi = 600,
       width = 8, height = 2.25)

#### loading in austin data: rr ####

load(file = "data/rpf.RData")

aus_rr = rpf %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 67329),
         r_wht_rate = (r_wht / 453034),
         r_lat_rate = (r_lat / 312603)) %>% 
  mutate(rr_bw = r_blk_rate / r_wht_rate,
         rr_lw = r_lat_rate / r_wht_rate) 

aus_rr$blmrv = 
  aus_rr$date - min(aus_rr$date[aus_rr$blm == 1])

#### loading in los angeles data: rr ####

load("data/stopdata_la.RData")

# rate ratio (pedestrian)

la_rr_ped = stopdat1 %>% 
  filter(reason_for_stop != "Traffic violation") %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 389998),
         r_wht_rate = (r_wht / 1169995),
         r_lat_rate = (r_lat / 1890299)) %>% 
  mutate(ped_rr_bw = r_blk_rate / r_wht_rate,
         ped_rr_lw = r_lat_rate / r_wht_rate) %>% 
  mutate(ped_rr_bw = ifelse(ped_rr_bw == Inf | ped_rr_bw == -Inf, NA, ped_rr_bw),
         ped_rr_lw = ifelse(ped_rr_lw == Inf | ped_rr_lw == -Inf, NA, ped_rr_lw))

la_rr_ped$blmrv = 
  la_rr_ped$date - min(la_rr_ped$date[la_rr_ped$blm == 1])

# rate ratio (vehicle)

la_rr_veh = stopdat1 %>% 
  filter(reason_for_stop == "Traffic violation") %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 389998),
         r_wht_rate = (r_wht / 1169995),
         r_lat_rate = (r_lat / 1890299)) %>% 
  mutate(veh_rr_bw = r_blk_rate / r_wht_rate,
         veh_rr_lw = r_lat_rate / r_wht_rate) %>% 
  mutate(veh_rr_bw = ifelse(veh_rr_bw == Inf | veh_rr_bw == -Inf, NA, veh_rr_bw),
         veh_rr_lw = ifelse(veh_rr_lw == Inf | veh_rr_lw == -Inf, NA, veh_rr_lw))

la_rr_veh$blmrv = 
  la_rr_veh$date - min(la_rr_veh$date[la_rr_veh$blm == 1])

#### loading in phl data: rr ####

# loading in data

load(file = "data/ph_stops_ped.RData")
load(file = "data/ph_stops_veh.RData")

# rate ratio (pedestrian)

phl_rr_ped = stops_ped %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 614254),
         r_wht_rate = (r_wht / 550102),
         r_lat_rate = (r_lat / 238965)) %>% 
  mutate(ped_rr_bw = r_blk_rate / r_wht_rate,
         ped_rr_lw = r_lat_rate / r_wht_rate) %>% 
  mutate(ped_rr_bw = ifelse(ped_rr_bw == Inf | ped_rr_bw == -Inf, NA, ped_rr_bw),
         ped_rr_lw = ifelse(ped_rr_lw == Inf | ped_rr_lw == -Inf, NA, ped_rr_lw))

phl_rr_ped$blmrv = 
  phl_rr_ped$date - min(phl_rr_ped$date[phl_rr_ped$blm == 1])

# rate ratio (vehicle)

phl_rr_veh = stops_veh %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 614254),
         r_wht_rate = (r_wht / 550102),
         r_lat_rate = (r_lat / 238965)) %>% 
  mutate(veh_rr_bw = r_blk_rate / r_wht_rate,
         veh_rr_lw = r_lat_rate / r_wht_rate) %>% 
  mutate(veh_rr_bw = ifelse(veh_rr_bw == Inf | veh_rr_bw == -Inf, NA, veh_rr_bw),
         veh_rr_lw = ifelse(veh_rr_lw == Inf | veh_rr_lw == -Inf, NA, veh_rr_lw))

phl_rr_veh$blmrv = 
  phl_rr_veh$date - min(phl_rr_veh$date[phl_rr_veh$blm == 1])

#### loading in sea data: rr ####

load("data/terry.RData")

terry = terry %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE", 1, 0))

sea_rr = 
  terry %>% 
  mutate(count = 1,
         r_unk = ifelse(subject_race == "Unknown" | subject_race == "-", 1, 0),
         r_wht = ifelse(subject_race == "White", 1, 0),
         r_asn = ifelse(subject_race == "Asian", 1, 0),
         r_blk = ifelse(subject_race == "Black or African American", 1, 0),
         r_lat = ifelse(subject_race == "Hispanic", 1, 0)) %>% 
  filter(initiate_off == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            arrest = mean(arrest, na.rm = TRUE),
            r_lat = sum(r_lat, na.rm = TRUE),
            r_wht = sum(r_wht, na.rm = TRUE),
            r_blk = sum(r_blk, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2019-06-01")) %>% 
  mutate(wht_rate = (r_wht / (.595 * 737015)) * 10000,
         blk_rate = (r_blk / (.068 * 737015)) * 10000,
         lat_rate = (r_lat / (.082 * 737015)) * 10000) %>% 
  mutate(rr_bw = blk_rate / wht_rate) %>% 
  mutate(rr_lw = lat_rate / wht_rate) %>%
  mutate(rr_bw = ifelse(rr_bw == Inf | is.na(rr_bw), NA, rr_bw)) %>% 
  mutate(rr_lw = ifelse(rr_lw == Inf | is.na(rr_lw), NA, rr_lw)) 

sea_rr$blmrv = 
  sea_rr$date - min(sea_rr$date[sea_rr$blm == 1])

#### analysis: rr ####

# austin

m_aus_rr_traffic_bw = 
  with(aus_rr, rdrobust(x = blmrv, y = rr_bw,
                        kernel = "uni", p = 1, all = TRUE))

m_aus_rr_traffic_lw = 
  with(aus_rr, rdrobust(x = blmrv, y = rr_lw,
                        kernel = "uni", p = 1, all = TRUE))

# los angeles

m_la_rr_ped_bw = 
  with(la_rr_ped, rdrobust(x = blmrv, y = ped_rr_bw,
                           kernel = "uni", p = 1, all = TRUE))

m_la_rr_ped_lw = 
  with(la_rr_ped, rdrobust(x = blmrv, y = ped_rr_lw,
                           kernel = "uni", p = 1, all = TRUE))

m_la_rr_traffic_bw = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_bw,
                           kernel = "uni", p = 1, all = TRUE))

m_la_rr_traffic_lw = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "uni", p = 1, all = TRUE))

# philadelphia

m_phl_rr_ped_bw = 
  with(phl_rr_ped, rdrobust(x = blmrv, y = ped_rr_bw,
                            kernel = "uni", p = 1, all = TRUE))

m_phl_rr_ped_lw = 
  with(phl_rr_ped, rdrobust(x = blmrv, y = ped_rr_lw,
                            kernel = "uni", p = 1, all = TRUE))

m_phl_rr_traffic_bw = 
  with(phl_rr_veh, rdrobust(x = blmrv, y = veh_rr_bw,
                            kernel = "uni", p = 1, all = TRUE))

m_phl_rr_traffic_lw = 
  with(phl_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                            kernel = "uni", p = 1, all = TRUE))

# seattle

m_sea_rr_terry_bw = 
  with(sea_rr, rdrobust(x = blmrv, y = rr_bw,
                        kernel = "uni", p = 1, all = TRUE))

#m_sea_rr_terry_lw = 
#  with(sea_rr, rdrobust(x = blmrv, y = rr_lw,
 #                       kernel = "uni", p = 1, all = TRUE))


df_rr_est = data.frame(
  
  est = c(m_aus_rr_traffic_bw$coef[1],
          m_aus_rr_traffic_lw$coef[1],
          m_la_rr_ped_bw$coef[1],
          m_la_rr_ped_lw$coef[1],
          m_la_rr_traffic_bw$coef[1],
          m_la_rr_traffic_lw$coef[1],
          m_phl_rr_ped_bw$coef[1],
          m_phl_rr_ped_lw$coef[1],
          m_phl_rr_traffic_bw$coef[1],
          m_phl_rr_traffic_lw$coef[1],
          m_sea_rr_terry_bw$coef[1]),
  
  
  se = c(m_aus_rr_traffic_bw$se[3],
         m_aus_rr_traffic_lw$se[3],
         m_la_rr_ped_bw$se[3],
         m_la_rr_ped_lw$se[3],
         m_la_rr_traffic_bw$se[3],
         m_la_rr_traffic_lw$se[3],
         m_phl_rr_ped_bw$se[3],
         m_phl_rr_ped_lw$se[3],
         m_phl_rr_traffic_bw$se[3],
         m_phl_rr_traffic_lw$se[3],
         m_sea_rr_terry_bw$se[3]),
  
  pv = c(m_aus_rr_traffic_bw$pv[3],
         m_aus_rr_traffic_lw$pv[3],
         m_la_rr_ped_bw$pv[3],
         m_la_rr_ped_lw$pv[3],
         m_la_rr_traffic_bw$pv[3],
         m_la_rr_traffic_lw$pv[3],
         m_phl_rr_ped_bw$pv[3],
         m_phl_rr_ped_lw$pv[3],
         m_phl_rr_traffic_bw$pv[3],
         m_phl_rr_traffic_lw$pv[3],
         m_sea_rr_terry_bw$pv[3]),
  
  city = c("AUS", "AUS", rep("LA", 4), rep("PHL", 4),
           "SEA"),
  
  outcome = factor(c("B/W RR (Traffic)", "L/W RR (Traffic)",
                     "B/W RR (Pedestrian)", "L/W RR (Pedestrian)",
                     "B/W RR (Traffic)", "L/W RR (Traffic)",
                     "B/W RR (Pedestrian)", "L/W RR (Pedestrian)",
                     "B/W RR (Traffic)", "L/W RR (Traffic)",
                     "B/W RR (Terry)"),
                   levels = c("B/W RR (Traffic)", "L/W RR (Traffic)",
                              "B/W RR (Pedestrian)", "L/W RR (Pedestrian)",
                              "B/W RR (Terry)"))
  
)

df_rr_est_plot = 
  df_rr_est %>% 
  filter(!grepl(x = outcome, "estrian")) 

meta_res_bw = 
  metagen(TE = (df_rr_est_plot %>% filter(grepl(x = outcome, "^B")))$est,
          seTE = (df_rr_est_plot %>% filter(grepl(x = outcome, "^B")))$se, 
          sm = "MD",
          method.random.ci = 'HK')

meta_res_lw = 
  metagen(TE = (df_rr_est_plot %>% filter(grepl(x = outcome, "^L")))$est,
          seTE = (df_rr_est_plot %>% filter(grepl(x = outcome, "^L")))$se, 
          sm = "MD",
          method.random.ci = 'HK')


df_hit_rate_plot = 
  bind_rows(df_rr_est_plot, data.frame(
    est = meta_res_bw$TE.random,
    se = meta_res_bw$seTE.random,
    pv = meta_res_bw$pval.random,
    city = "B/W Meta-\nAnalysis",
    outcome = "B/W Meta-Analysis"
  ),
  data.frame(
    est = meta_res_lw$TE.random,
    se = meta_res_lw$seTE.random,
    pv = meta_res_lw$pval.random,
    city = "L/W Meta-\nAnalysis",
    outcome = "L/W Meta-Analysis"
  )) %>% 
  mutate(city = factor(city,
                       levels = c("AUS", "LA", "PHL", "SEA",
                                  "B/W Meta-\nAnalysis",
                                  "L/W Meta-\nAnalysis")),
         outcome = factor(outcome,
                          levels = c("B/W RR (Traffic)",
                                     "L/W RR (Traffic)",
                                     "B/W RR (Terry)",
                                     "B/W Meta-Analysis",
                                     "L/W Meta-Analysis")))

segment3 = 
  df_hit_rate_plot %>% 
  ggplot() + 
  geom_hline(yintercept = 0,
             size = .2) + 
  geom_point(aes(x = city,
                 y = est,
                 shape = outcome),
             position = position_dodge(.25),
             size = 1.25) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                position = position_dodge(.25),
                width = 0,
                size = .4) +
  labs(x = "City", y = "Coefficient\n(BLM Protest)",
       title = "C. Rate Ratios",
       shape = "Outcome") + 
  theme_tufte(base_size = 9) +
  theme(legend.key.size = unit(.1, 'cm'),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 4),
        axis.text.x = element_text(size = 5),
        legend.position = "bottom") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE))

grob_segments = arrangeGrob(segment1, segment2, segment3, ncol = 3)

ggsave(plot = grob_segments, filename = "pics/figure3res.jpeg", 
       width = 8, height = 2.25,
       dpi = 600)

# testing robustness of LA finding 

m_la_rr_traffic_lw1 = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "uni", p = 1, all = TRUE))
m_la_rr_traffic_lw2 = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "uni", p = 2, all = TRUE))
m_la_rr_traffic_lw3 = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "uni", p = 3, all = TRUE))

m_la_rr_traffic_lw4 = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "tri", p = 1, all = TRUE))
m_la_rr_traffic_lw5 = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "tri", p = 2, all = TRUE))
m_la_rr_traffic_lw6 = 
  with(la_rr_veh, rdrobust(x = blmrv, y = veh_rr_lw,
                           kernel = "tri", p = 3, all = TRUE))

# pretty robust 

df_output_la = matrix(NA, nrow = 100, ncol = 4) %>% 
  as.data.frame %>% 
  `colnames<-` (c("est", 'se', 'pv', 'dayscut'))

for (i in 1:100) {
  
  print(paste0("Iteration ", i))
  
  la_rr_veh_cutd = 
    (la_rr_veh %>% filter(blmrv >= 1 | blmrv <= -1) %>% 
       mutate(blmrv = ifelse(blmrv > 0, blmrv - i, blmrv)))
  
  outremove = 
    with(la_rr_veh_cutd, rdrobust(x = blmrv, y = veh_rr_lw,
                                kernel = "uni", p = 1, all = TRUE))
  
  df_output_la$est[i] = outremove$coef[1]
  df_output_la$se[i] = outremove$se[3]
  df_output_la$pv[i] = outremove$pv[3]
  df_output_la$dayscut[i] = i
  
}

df_output_la %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = dayscut, y = est)) + 
  geom_errorbar(aes(x = dayscut, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                width = 0) + 
  labs(x = "Days Cut", y = "RDiT Coefficient") + 
  theme_tufte()

#### placebos before and after: black/white ratio #### 

phl_rr_veh$rr_bw = phl_rr_veh$veh_rr_bw
list_to_loop_df = 
  list(phl_rr_veh, sea_rr)
city_label = c("PHL", "SEA")
list_to_loop_df_out = as.list(rep(NA, length(list_to_loop_df)))
seq_2_l = seq(from = -100, to = 100, by = 20)

for (i in 1:length(list_to_loop_df)) {
  
  print(paste0("I iteration ", i))
  
  df_seq_out = 
    matrix(NA, nrow = length(seq_2_l), ncol = 6) %>% 
    as.data.frame %>% 
    `colnames<-` (c("est", "se", "pv", 'plcval', 'city', 'outcome'))
  
  for (k in 1:length(seq_2_l)) {
    
    print(paste0("K iteration ", k))
    
    out = with(list_to_loop_df[[i]],
               rdrobust(x = I(blmrv + seq_2_l[k]), y = rr_bw, p = 1,
                        kernel = 'uni'))
    
    df_seq_out[k, 1] = out$coef[3]
    df_seq_out[k, 2] = out$se[3]
    df_seq_out[k, 3] = out$pv[3]
    df_seq_out[k, 4] = seq_2_l[k]
    df_seq_out[k, 5] = city_label[i]
    df_seq_out[k, 6] = stop_label[i]
    
  }
  
  list_to_loop_df_out[[i]] = df_seq_out
  
}

list_to_loop_df_out %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(city_stop = paste0(city, ', ', outcome)) %>% 
  mutate(plcvalbin = ifelse(plcval == 0, 1, 0)) %>% 
  mutate(plcvalbin = factor(plcvalbin, levels = c(0, 1))) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = plcval, y = est,
                 col = plcvalbin),
             size = 2) + 
  geom_errorbar(aes(x = plcval, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = plcvalbin),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = plcval, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = plcvalbin),
                size = .8,
                width = 0) + 
  guides(col = 'none') + 
  labs(x = "Placebo Discontinuity (Days to BLM Protest)",
       y = "RDiT Coefficient",
       title = 'Placebo Discontinuities Before and After the BLM Protest (Black/White Rate Ratio Outcome)') + 
  scale_color_grey(start = .6, end = 0) + 
  facet_wrap(~city_stop, scales = 'free') + 
  theme_tufte()

ggsave(plot = last_plot(), filename = 'pics/plcdisc_rratio.jpeg', dpi = 600,
       width = 8, height = 2.5)

#### MAKE REGRESSION TABLE FOR APPENDIX ####

library(xtable)

# hit rate rows
atx_hit <- as.data.frame(cbind(m_aus_hit_traffic$coef[1],m_aus_hit_traffic$se[3],
                               m_aus_hit_traffic$pv[3],m_aus_hit_traffic$N[1],
                               (m_aus_hit_traffic$bws[1]*2),m_aus_hit_traffic$bws[1]))

la_hit <- as.data.frame(cbind(m_la_hit_traffic$coef[1],m_la_hit_traffic$se[3],
                              m_la_hit_traffic$pv[3],m_la_hit_traffic$N[1],
                              (m_la_hit_traffic$bws[1]*2),m_la_hit_traffic$bws[1]))

phl_hit <- as.data.frame(cbind(m_phl_hit_traffic$coef[1],m_phl_hit_traffic$se[3],
                               m_phl_hit_traffic$pv[3],m_phl_hit_traffic$N[1],
                               (m_phl_hit_traffic$bws[1]*2),m_phl_hit_traffic$bws[1]))

sea_hit <- as.data.frame(cbind(m_sea_hit$coef[1],m_sea_hit$se[3],
                               m_sea_hit$pv[3],m_sea_hit$N[1],
                               (m_sea_hit$bws[1]*2),m_sea_hit$bws[1]))

# arrest rate rows
atx_arrest <- as.data.frame(cbind(m_aus_arr_traffic$coef[1],m_aus_arr_traffic$se[3],
                                  m_aus_arr_traffic$pv[3],m_aus_arr_traffic$N[1],
                                  (m_aus_arr_traffic$bws[1]*2),m_aus_arr_traffic$bws[1]))

la_arrest <- as.data.frame(cbind(m_la_arr_traffic$coef[1],m_la_arr_traffic$se[3],
                                 m_la_arr_traffic$pv[3],m_la_arr_traffic$N[1],
                                 (m_la_arr_traffic$bws[1]*2),m_la_arr_traffic$bws[1]))

phl_arrest <- as.data.frame(cbind(m_phl_arr_traffic$coef[1],m_phl_arr_traffic$se[3],
                                  m_phl_arr_traffic$pv[3],m_phl_arr_traffic$N[1],
                                  (m_phl_arr_traffic$bws[1]*2),m_phl_arr_traffic$bws[1]))

sea_arrest <- as.data.frame(cbind(m_sea_arr$coef[1],m_sea_arr$se[3],
                                  m_sea_arr$pv[3],m_sea_arr$N[1],
                                  (m_sea_arr$bws[1]*2),m_sea_arr$bws[1]))

# b/w rate ratio rows

atx_ratio <- as.data.frame(cbind(m_aus_rr_traffic_bw$coef[1],m_aus_rr_traffic_bw$se[3],
                                 m_aus_rr_traffic_bw$pv[3],m_aus_rr_traffic_bw$N[1],
                                 (m_aus_rr_traffic_bw$bws[1]*2),m_aus_rr_traffic_bw$bws[1]))

la_ratio <- as.data.frame(cbind(m_la_rr_traffic_bw$coef[1],m_la_rr_traffic_bw$se[3],
                                m_la_rr_traffic_bw$pv[3],m_la_rr_traffic_bw$N[1],
                                (m_la_rr_traffic_bw$bws[1]*2),m_la_rr_traffic_bw$bws[1]))

phl_ratio <- as.data.frame(cbind(m_phl_rr_traffic_bw$coef[1],m_phl_rr_traffic_bw$se[3],
                                 m_phl_rr_traffic_bw$pv[3],m_phl_rr_traffic_bw$N[1],
                                 (m_phl_rr_traffic_bw$bws[1]*2),m_phl_rr_traffic_bw$bws[1]))

sea_ratio <- as.data.frame(cbind(m_sea_rr_terry_bw$coef[1],m_sea_rr_terry_bw$se[3],
                                 m_sea_rr_terry_bw$pv[3],m_sea_rr_terry_bw$N[1],
                                 (m_sea_rr_terry_bw$bws[1]*2),m_sea_rr_terry_bw$bws[1]))


tab <- as.data.frame(rbind(atx_hit,la_hit,phl_hit,sea_hit,
                           atx_arrest,la_arrest,phl_hit,sea_hit,
                           atx_ratio,la_ratio,phl_ratio,sea_ratio))

tab <- tab %>% mutate(V4 = as.character(V4))
tab <- tab %>% mutate_if(is.numeric,round, digits = 2)

cities <- as.data.frame(rbind("Austin","LA","Philly","Seattle",
                              "Austin","LA","Philly","Seattle",
                              "Austin","LA","Philly","Seattle"))
outcome <- as.data.frame(rbind("Hit rate","Hit rate","Hit rate","Hit rate",
                               "Arrest rate","Arrest rate","Arrest rate","Arrest rate",
                               "B/W rate ratio","B/W rate ratio","B/W rate ratio","B/W rate ratio"))

tab <- as.data.frame(cbind(cities,outcome,tab))
colnames(tab) <- c("City","Outcome","Coeff","SE","P-Val","N-Val","Effective N","Bandwidth Est.")

xt <- xtable(
  x = tab,
  align = "lllcccccc",
  type = "html",
  caption = "RDiT Coefficients Characterizing the Effect of BLM Protests on Policing Quality.",
  label = "tab:quality_reg_tab"
)


addtorow <- list()
addtorow$pos <- list(nrow(xt))
addtorow$command <- as.vector(c(paste("\\hline\n\\multicolumn{8}{l}{\\emph{All estimates are specified with a uniform kernel and polynomial degree equal to 1.}}\\\\\n",
                                      paste("\\multicolumn{8}{l}{\\emph{Standard errors are robust.}}\\\\\n",
                                            sep=""))))

print(xt,
      format.args = list(big.mark = ","),
      caption.placement = "top",
      include.rownames = FALSE,
      add.to.row = addtorow,
      file = "fig4_reg_tab.tex",
      hline.after = c(-1, 0, 12)
)


